/*****************************************************************************/
/**
 *  @file   Import.cpp
 *  @author Naohisa Sakamoto
 */
/*----------------------------------------------------------------------------
 *
 *  Copyright (c) Visualization Laboratory, Kyoto University.
 *  All rights reserved.
 *  See http://www.viz.media.kyoto-u.ac.jp/kvs/copyright/ for details.
 *
 *  $Id$
 */
/*****************************************************************************/
#include "Import.h"
#include "SPHFile.h"
#include "SPHFileList.h"
#include <kvs/AnyValueArray>
#include <kvs/ValueArray>
#include <kvs/Assert>

namespace
{

struct CombinedData
{
    kvs::Real32 min_value; ///< min. value
    kvs::Real32 max_value; ///< max. value
    kvs::AnyValueArray values; ///< value array
};

kvs::AnyValueArray ByteArray( const kvs::AnyValueArray& data )
{
    kvs::ValueArray<kvs::UInt8> values( data.size() );
    for ( size_t i = 0; i < values.size(); i++ )
    {
        const kvs::Real64 value = kvs::Math::Clamp( data.at<kvs::Real64>(i), 0.0, 1.0 );
        values[i] = kvs::Math::Round( 255.0 * value );
    }

    return values;
}

CombinedData CombinedByteArray( const kvsbloodvis::util::SPHFileList& sph_list )
{
    KVS_ASSERT( sph_list.svType() == 1 );

    const kvs::Vector3d global_org = sph_list.org();
    const kvs::Vector3d global_pitch = sph_list.pitch();
    const kvs::Vector3d global_max = sph_list.max();
    const kvs::Vector3ui resolution = sph_list.resolution();
    const kvs::UInt32 values_size = resolution.x() * resolution.y() *resolution.z();
    kvs::ValueArray<kvs::UInt8> values( values_size );
    for ( size_t index = 0; index < sph_list.size(); index++ )
    {
        const kvs::Vector3d local_max = sph_list.sph( index ).max();
        const kvs::Vector3d local_org = sph_list.sph( index ).org();
        const kvs::Vector3d local_pitch = sph_list.sph( index ).pitch();
        const kvs::AnyValueArray local_data = sph_list.sph( index ).data();
        const size_t I_start = ( local_org[0] - global_org[0] ) / global_pitch[0];
        const size_t J_start = ( local_org[1] - global_org[1] ) / global_pitch[1];
        const size_t K_start = ( local_org[2] - global_org[2] ) / global_pitch[2];
        const size_t i_end = local_max[0];
        const size_t j_end = local_max[1];
        const size_t k_end = local_max[2];
        for ( size_t k = 0, K = K_start; k < k_end; k++, K++ )
        {
            for ( size_t j = 0, J = J_start; j < j_end; j++, J++ )
            {
                for ( size_t i = 0, I = I_start; i < i_end; i++, I++ )
                {
                    const size_t local_index = k * local_max[0] * local_max[1] + j * local_max[0] + i;
                    const size_t global_index = K * global_max[0] * global_max[1] + J * global_max[0] + I;
                    kvs::Real64 value = kvs::Math::Clamp( local_data.at<kvs::Real64>( local_index ), 0.0, 1.0 );
                    values[ global_index ] = kvs::Math::Round( 255.0 * value );
                }
            }
        }
    }

    CombinedData data;
    data.min_value = kvs::Value<kvs::UInt8>::Min();
    data.max_value = kvs::Value<kvs::UInt8>::Max();
    data.values = values;

    return data;
}

CombinedData CombinedFloatArray( const kvsbloodvis::util::SPHFileList& sph_list )
{
    const size_t veclen = sph_list.svType() == 1 ? 1 : 3;
    const kvs::Vector3ui resolution = sph_list.resolution();

    const kvs::Vector3d global_org = sph_list.org();
    const kvs::Vector3d global_pitch = sph_list.pitch();
    const kvs::Vector3d global_max = sph_list.max();
    const kvs::UInt32 values_size = resolution.x() * resolution.y() *resolution.z();
    kvs::ValueArray<kvs::Real32> values( values_size * veclen );
    values.fill(0);

    kvs::Real32 min_value = kvs::Value<kvs::Real32>::Max();
    kvs::Real32 max_value = kvs::Value<kvs::Real32>::Min();

    for ( size_t index = 0; index < sph_list.size(); index++ )
    {
        const kvs::Vector3d local_max = sph_list.sph( index ).max();
        const kvs::Vector3d local_org = sph_list.sph( index ).org();
        const kvs::Vector3d local_pitch = sph_list.sph( index ).pitch();
        const kvs::AnyValueArray local_data = sph_list.sph( index ).data();
        const size_t I_start = ( local_org[0] - global_org[0] ) / global_pitch[0];
        const size_t J_start = ( local_org[1] - global_org[1] ) / global_pitch[1];
        const size_t K_start = ( local_org[2] - global_org[2] ) / global_pitch[2];
        const size_t i_end = local_max[0];
        const size_t j_end = local_max[1];
        const size_t k_end = local_max[2];
        for ( size_t k = 0, K = K_start; k < k_end; k++, K++ )
        {
            for ( size_t j = 0, J = J_start; j < j_end; j++, J++ )
            {
                for ( size_t i = 0, I = I_start; i < i_end; i++, I++ )
                {
                    const size_t local_index = k * local_max[0] * local_max[1] + j * local_max[0] + i;
                    const size_t global_index = K * global_max[0] * global_max[1] + J * global_max[0] + I;
                    kvs::Real32 temp_value = 0.0f;
                    for ( size_t v = 0; v < veclen; v++ )
                    {
                        const kvs::Real32 value = local_data.at<kvs::Real32>( veclen * local_index + v );
                        values[ veclen * global_index + v ] = value;
                        temp_value += value * value;
                    }

                    temp_value = std::sqrt( temp_value );
                    min_value = kvs::Math::Min( min_value, temp_value );
                    max_value = kvs::Math::Max( max_value, temp_value );
                }
            }
        }
    }

    CombinedData data;
    data.min_value = min_value;
    data.max_value = max_value;
    data.values = values;

    return data;
}

kvs::Real32 Extrap(
    const size_t index0,
    const size_t index1,
    const size_t stride,
    const kvs::AnyValueArray& data )
{
    const size_t veclen = 3;
    const kvs::Real32 s0 = data.at<kvs::Real32>( index0 * veclen + stride );
    const kvs::Real32 s1 = data.at<kvs::Real32>( index1 * veclen + stride );
    return 1.5f * s0 - 0.5f * s1;
}

kvs::Real32 Interp(
    const size_t index0,
    const size_t index1,
    const size_t stride,
    const kvs::AnyValueArray& data )
{
    const size_t veclen = 3;
    const kvs::Real32 s0 = data.at<kvs::Real32>( index0 * veclen + stride );
    const kvs::Real32 s1 = data.at<kvs::Real32>( index1 * veclen + stride );
    return 0.5f * ( s0 + s1 );
}

kvs::ValueArray<kvs::Real32> StaggeredToColocate( const kvs::AnyValueArray& data, const kvs::Vector3ui& dim )
{
    const size_t veclen = 3;
    const size_t value_size = dim.x() * dim.y() * dim.z() * veclen;
    kvs::ValueArray<kvs::Real32> values( value_size );

    const size_t nnodes_per_slice = dim.x() * dim.y();
    const size_t nnodes_per_line = dim.x();
    for ( size_t k = 0, index = 0; k < dim.z(); k++ )
    {
        for ( size_t j = 0; j < dim.y(); j++ )
        {
            for ( size_t i = 0; i < dim.x(); i++, index++ )
            {
                const size_t u0 = i == 0 ? index : index - 1;
                const size_t u1 = i == 0 ? index + 1 : index;
                const kvs::Real32 u = i == 0 ? Extrap( u0, u1, 0, data ) : Interp( u0, u1, 0, data );

                const size_t v0 = j == 0 ? index : index - nnodes_per_line;
                const size_t v1 = j == 0 ? index + nnodes_per_line : index;
                const kvs::Real32 v = j == 0 ? Extrap( v0, v1, 1, data ) : Interp( v0, v1, 1, data );

                const size_t w0 = k == 0 ? index : index - nnodes_per_slice;
                const size_t w1 = k == 0 ? index + nnodes_per_slice : index;
                const kvs::Real32 w = k == 0 ? Extrap( w0, w1, 2, data ) : Interp( w0, w1, 2, data );

                values[ index * veclen + 0 ] = u;
                values[ index * veclen + 1 ] = v;
                values[ index * veclen + 2 ] = w;
            }
        }
    }

    return values;
}

}

namespace kvsbloodvis
{

namespace util
{

kvs::StructuredVolumeObject Import( const SPHFile& sph, const bool byte )
{
    const size_t veclen = sph.svType() == 1 ? 1 : 3;
    const kvs::Vector3ui resolution = sph.resolution();

    kvs::AnyValueArray values;
    if ( veclen == 1 )
    {
        if ( byte ) values = ::ByteArray( sph.data() );
        else values = sph.data();
    }
    else // veclen == 3
    {
        values = ::StaggeredToColocate( sph.data(), resolution );
    }

    kvs::StructuredVolumeObject volume;
    volume.setGridType( kvs::StructuredVolumeObject::Uniform );
    volume.setVeclen( veclen );
    volume.setResolution( resolution );
    volume.setValues( values );
    volume.updateMinMaxCoords();
    volume.updateMinMaxValues();

    return volume;
}

kvs::StructuredVolumeObject Import( const SPHFileList& sph_list, const bool byte )
{
    const size_t veclen = sph_list.svType() == 1 ? 1 : 3;
    const kvs::Vector3ui resolution = sph_list.resolution();

    CombinedData data;
    if ( veclen == 1 )
    {
        if ( byte ) data = ::CombinedByteArray( sph_list );
        else data = ::CombinedFloatArray( sph_list );
    }
    else // veclen == 3
    {
        data = ::CombinedFloatArray( sph_list );
        data.values = ::StaggeredToColocate( data.values, resolution );
    }

    kvs::StructuredVolumeObject volume;
    volume.setGridType( kvs::StructuredVolumeObject::Uniform );
    volume.setVeclen( veclen );
    volume.setResolution( resolution );
    volume.setValues( data.values );
    volume.setMinMaxValues( data.min_value, data.max_value );
    volume.updateMinMaxCoords();

    return volume;
}

} // end of namespace util

} // end of namespace kvsbloodvis
